Nonlinear Regression (Archontoulis and Miguez) paper

https://femiguez.github.io/nlraa-docs/nlraa-Oddi-LFMC.html

https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/ece3.5543

#install.packages("emmeans")

library(nlme)
library(nlraa)
library(ggplot2)
library(emmeans)   # emmeans::contrast(), emmeans::emmeans()

Steps in fitting a nonlinear mixed model

Part I.

data(lfmc)
lfmc.gd <- groupedData(lfmc ~ time | group, data = lfmc, order.groups = FALSE) 
lfmc.gd$plot = with(lfmc.gd, factor(plot, levels =c("4", "5", "6", "1", "2", "3")))
head(lfmc.gd)

Grouped Data: lfmc ~ time | group
         leaf.type time plot site  lfmc     group
1      M. spinosum    1    1    P 238.4 SM plot 1
2      M. spinosum    1    1    P 269.8 SM plot 1
3      M. spinosum    1    1    P 325.2 SM plot 1
4 S. bracteolactus    1    1    P 326.7 SS plot 1
5 S. bracteolactus    1    1    P 257.3 SS plot 1
6 S. bracteolactus    1    1    P 279.3 SS plot 1

ggplot(lfmc, aes(time, lfmc)) +
  geom_point() +
  geom_smooth(method = "loess") +
  facet_wrap(~leaf.type) +
  xlab("Time (days)") +
  ylab("LMFC (%)")

`geom_smooth()` using formula 'y ~ x'

NONLINEAR MIXED-EFFECTS MODEL (M1)

nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)

Warning: 3 errors caught in nls(model, data = data, control = controlvals).  The error messages and their frequencies are

step factor 0.000488281 reduced below 'minFactor' of 0.000976562 
                                                               1 
                     number of iterations exceeded maximum of 50 
                                                               2 

coef(nlsL)

plot(intervals(nlsL))

Random-Effects Mode

nl.re <- nlme(nlsL, control = nlmeControl(maxIter = 5000, msMaxIter = 1500))

Warning in (function (model, data = sys.frame(sys.parent()), fixed, random,  :
  Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: function evaluation limit reached without convergence (9)

nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
fxf <- fixef(nl.re.1)
fxf

        A         w         m         s 
193.45254  21.45745  31.48273 -22.59708 

# nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
#                 start = c(A = fxf[1],0,0,0,
#                           w = fxf[2],0,0,0,
#                           m = fxf[3],0,0,0,
#                           s = fxf[4],0,0,0))

nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
fxf <- fixef(nl.re.2)
fxf

        A         w         m         s 
193.46654  21.43199  31.48852 -22.60602 

nl.me.1 <- update(
  nl.re.2, 
  fixed = list(A + w + m + s ~ leaf.type),  
  start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4],0,0,0))
nl.me.1

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1034.375
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E 
                  23.135414                   53.294908 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                 296.620596                  263.313118 
              w.(Intercept)          w.leaf.typeGrass E 
                  50.291577                  -34.313924 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                 -26.732287                    2.135665 
              m.(Intercept)          m.leaf.typeGrass E 
                  51.533983                  -21.920404 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                 -20.647403                  -18.170359 
              s.(Intercept)          s.leaf.typeGrass E 
                  17.605794                  -25.926031 
     s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                 -43.529450                  -33.666226 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
        A.(Intercept) w.(Intercept) m.(Intercept) Residual
StdDev:      8.112813   0.001245815      2.468334 15.38607

Number of Observations: 247
Number of Groups: 12 

nl.re.2

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1070.648
  Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1) 
        A         w         m         s 
193.46654  21.43199  31.48852 -22.60602 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
               A        w        m Residual
StdDev: 119.3792 10.94165 8.927279 15.26574

Number of Observations: 247
Number of Groups: 12 

AIC(nl.re.2, nl.me.1)

IC_tab(nl.re.2, nl.me.1)

nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
IC_tab(nl.me.1, nl.me.2)

nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1)) 
IC_tab(nl.me.2, nl.me.3)

plot(nl.me.3)

rse <- residuals(nl.me.3) / sigma(nl.me.3)
hist(residuals(nl.me.3, type = "normalized"))

#par(pty = "s")
qqnorm(rse)
qqline(rse)

nl.me.3.vm <- update(nl.me.3, weights = varIdent(form = ~1 | leaf.type))
IC_tab(nl.me.3, nl.me.3.vm)

plot(nl.me.3.vm)

hist(residuals(nl.me.3.vm, type = "normalized"), main = "", xlab = "Residuals")

qqnorm(resid(nl.me.3.vm, type = "normalized"))
qqline(resid(nl.me.3.vm, type = "normalized"))

plot(ACF(nl.me.3.vm), alpha = 0.01)

Evaluation of Fixed Effects

nl.me.ml <- update(nl.me.3.vm, method = "ML")
nl.me.ml

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -965.5628
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E 
                 50.1522850                  28.9553860 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                346.5816365                 231.1569668 
              w.(Intercept)          w.leaf.typeGrass E 
                 24.6693608                  -9.4837876 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                -35.8327959                  29.7533014 
              m.(Intercept)          m.leaf.typeGrass E 
                 48.9613646                 -19.9780200 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                -26.2920797                 -14.9909275 
              s.(Intercept)          s.leaf.typeGrass E 
                -16.2294126                   6.3520529 
     s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                -21.0664461                   0.9389322 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.412262 6.824811

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8547836        3.1456077        3.0262123 
Number of Observations: 247
Number of Groups: 12 

nl.me.ml.1 <- update(nl.me.ml,
  fixed = list(A + w + m ~ leaf.type, s ~ 1),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4])
)
nl.me.ml.1

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -969.2617
  Fixed: list(A + w + m ~ leaf.type, s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   49.94562                    39.60979 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  219.45652                   232.11909 
              w.(Intercept)          w.leaf.typeGrass E 
                   25.24501                   -14.34293 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   35.36121                    28.70388 
              m.(Intercept)          m.leaf.typeGrass E 
                   48.30390                   -21.59477 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                  -15.57765                   -14.37454 
                          s 
                  -15.44291 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.612472 6.825604

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8714808        3.2750695        3.0238813 
Number of Observations: 247
Number of Groups: 12 

IC_tab(nl.me.ml, nl.me.ml.1) 

nl.me.ml.2 <- update(nl.me.ml, fixed = list(A + w ~ leaf.type, m + s ~ 1),  
                     start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3], s=fxf[4])) 
nl.me.ml.2

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -972.7465
  Fixed: list(A + w ~ leaf.type, m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   54.25730                    30.93396 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  223.27914                   240.31543 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.05803                   -20.68663 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   31.68646                    26.99003 
                          m                           s 
                   30.91655                   -16.15435 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.216113 7.028345

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8820413        3.1693967        2.9643305 
Number of Observations: 247
Number of Groups: 12 

nl.me.ml.3 <- update(nl.me.ml, fixed = list(A ~ leaf.type, w + m + s  ~ 1),  
                     start = c(A=fxf[1],0,0,0, w=fxf[2], m=fxf[3], s=fxf[4])) 
nl.me.ml.3

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1017.91
  Fixed: list(A ~ leaf.type, w + m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   63.92214                    14.27256 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  251.99464                   265.42794 
                          w                           m 
                   17.33007                    32.65667 
                          s 
                  -25.16304 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:       8.13893 7.795662

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
        1.000000         1.584893         2.813972         2.731419 
Number of Observations: 247
Number of Groups: 12 

nl.me.ml.4 <- update(nl.me.ml, fixed = list(w ~ leaf.type, A + m + s  ~ 1),  
                     start = c(A=fxf[1], w=fxf[2],0,0,0, m=fxf[3], s=fxf[4])) 
nl.me.ml.4

Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -999.8224
  Fixed: list(w ~ leaf.type, A + m + s ~ 1) 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.16028                   -20.38030 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   33.19422                    28.72473 
                          A                           m 
                  175.78134                    31.08599 
                          s 
                  -15.70300 

Random effects:
 Formula: A ~ 1 | group
               A Residual
StdDev: 108.0655 7.042353

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8770644        3.1888994        2.9304282 
Number of Observations: 247
Number of Groups: 12 

IC_tab(nl.me.ml.1, nl.me.ml.2, nl.me.ml.3, nl.me.ml.4)

M1 <- update(nl.me.ml.2, method = "REML")
summary(M1)

Nonlinear mixed-effects model fit by REML
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.408521 7.162899

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8851762        3.1796110        2.9647390 
Fixed effects:  list(A + w ~ leaf.type, m + s ~ 1) 

 Correlation: 
                            A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE w..M.s
A.leaf.typeGrass E          -0.527                                          
A.leaf.typeM. spinosum      -0.146  0.535                                   
A.leaf.typeS. bracteolactus -0.115  0.537  0.746                            
w.(Intercept)               -0.217 -0.034 -0.200 -0.216                     
w.leaf.typeGrass E          -0.035 -0.283 -0.382 -0.399 -0.258              
w.leaf.typeM. spinosum      -0.118 -0.227 -0.547 -0.454  0.157  0.573       
w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575  0.188  0.601  0.640
m                           -0.217 -0.324 -0.633 -0.666  0.092  0.145  0.161
s                           -0.246 -0.354 -0.710 -0.748  0.416  0.575  0.702
                            w..S.b m     
A.leaf.typeGrass E                       
A.leaf.typeM. spinosum                   
A.leaf.typeS. bracteolactus              
w.(Intercept)                            
w.leaf.typeGrass E                       
w.leaf.typeM. spinosum                   
w.leaf.typeS. bracteolactus              
m                            0.172       
s                            0.752  0.544

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.27999716 -0.59657157 -0.02042043  0.57273326  3.12017921 

Number of Observations: 247
Number of Groups: 12 

fixef(M1)

              A.(Intercept)          A.leaf.typeGrass E 
                   54.25350                    30.92293 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  223.24504                   240.27654 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.05581                   -20.68938 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   31.67297                    26.97508 
                          m                           s 
                   30.92881                   -16.15406 

ranef(M1)

          A.(Intercept)
GW plot 4     -1.915443
GW plot 5     -1.923198
GW plot 6      3.838641
GE plot 1     15.652459
SM plot 1      6.967303
SS plot 1      9.166925
GE plot 2    -11.066197
SM plot 2     -5.358910
SS plot 2      2.442781
GE plot 3     -4.586263
SM plot 3     -1.608393
SS plot 3    -11.609707

intervals(M1)

Approximate 95% confidence intervals

 Fixed effects:
                                lower      est.     upper
A.(Intercept)                42.45370  54.25350  66.05331
A.leaf.typeGrass E           13.51301  30.92293  48.33286
A.leaf.typeM. spinosum      191.88318 223.24504 254.60691
A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
w.(Intercept)                25.68039  29.05581  32.43122
w.leaf.typeGrass E          -25.79788 -20.68938 -15.58088
w.leaf.typeM. spinosum       16.39521  31.67297  46.95073
w.leaf.typeS. bracteolactus  11.07083  26.97508  42.87934
m                            26.85931  30.92881  34.99831
s                           -20.81510 -16.15406 -11.49302

 Random Effects:
  Level: group 
                     lower     est.    upper
sd(A.(Intercept)) 5.035417 9.408521 17.57953

 Variance function:
                     lower      est.    upper
Grass E          0.6772617 0.8851762 1.156919
M. spinosum      2.4528838 3.1796110 4.121649
S. bracteolactus 2.2864085 2.9647390 3.844316

 Within-group standard error:
   lower     est.    upper 
5.970704 7.162899 8.593144 

Part II

# 2.2.1 - Overall predictions (Table 3) ---- 
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW) 
Age <- fixef(M1)[1]+fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1]+fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1]+fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW) 
Wge <- fixef(M1)[5]+fixef(M1)[6] # "w" for grasses in the E site (GE) 
Wsm <- fixef(M1)[5]+fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5]+fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types

GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")

par(mfrow=c(2,2), cex=0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=2, lty=1, add=T, col="darkgreen") 

# Grasses E
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age-Wge)/(1 + exp((M-x)/S))+Wge), lwd=2, lty=1, add=T, col="green")

# M. Spinosum
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=2, lty=1, add=T, col="darkorange") 

# S. filaginoides 
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass-Wss)/(1 + exp((M-x)/S))+Wss), lwd=2, lty=1, add=T, col="darkred") 

# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1]+ranef(M1)[1,1] # "A" for grasses in the plot 4 (GW) 
Agw.p5 <- fixef(M1)[1]+ranef(M1)[2,1] # "A" for grasses in the plot 5 (GW) 
Agw.p6 <- fixef(M1)[1]+ranef(M1)[3,1] # "A" for grasses in the plot 6 (GW) 
Age.p1 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[4,1] # "A" for grasses in the plot 1 (GW) 
Age.p2 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[7,1] # "A" for grasses in the plot 2 (GW) 
Age.p3 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[10,1] # "A" for grasses in the plot 3 (GW) 
Asm.p1 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[5,1] # "A" for M. spinosum in the plot 1 (SM) 
Asm.p2 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[8,1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[11,1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[6,1] # "A" for S. filaginoides in the plot 1 (SM)   
Ass.p2 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[9,1] # "A" for S. filaginoides in the plot 2 (SM) 
Ass.p3 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[12,1] # "A" for S. filaginoides in the plot 3 (SM)
  
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
# Grasses W site 
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw.p4-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 4

curve(((Agw.p5-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 5
curve(((Agw.p6-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 6

# Grasses E site
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age.p1-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 1

curve(((Age.p2-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 2
curve(((Age.p3-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 3

# M. Spinosum (E site) 
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm.p1-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 1

curve(((Asm.p2-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 2
curve(((Asm.p3-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 3

# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass.p1-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 1

curve(((Ass.p2-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 2
curve(((Ass.p3-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 3

# 2.2.3 - Derivatives (drying speed) ---- 
lfmc.GW <- expression((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site

lfmc.GE <- expression((Age-Wge)/(1 + exp((M-x)/S))+Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site

lfmc.SM <- expression((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum

lfmc.SS <- expression((Ass-Wss)/(1 + exp((M-x)/S))+Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides

xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)

par(mfrow=c(2,2), cex=0.75) # (Figure 4)
plot(xvec, y.ds.GW, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the W site")
lines(xvec, -1*(attr(y.ds.GW, "grad")), lwd=2, lty=1, col="darkgreen")

plot(xvec, y.ds.GE, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the E site")
lines(xvec, -1*(attr(y.ds.GE, "grad")), lwd=2, lty=1, col="green")

plot(xvec, y.ds.SM, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="M. spinosum", font.main=4)
lines(xvec, -1*(attr(y.ds.SM, "grad")), lwd=2, lty=1, col="darkorange")

plot(xvec, y.ds.SS, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="S. filaginoides", font.main=4)
lines(xvec, -1*(attr(y.ds.SS, "grad")), lwd=2, lty=1, col="darkred")

# 2.3 ---- Testing differences among leaf types ------------------------------------------- 

# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")

 contrast                       estimate    SE  df t.ratio p.value
 Grass W - Grass E                 -30.9  8.84 226  -3.500  0.0031
 Grass W - M. spinosum            -223.2 15.92 226 -14.027  <.0001
 Grass W - S. bracteolactus       -240.3 16.86 226 -14.254  <.0001
 Grass E - M. spinosum            -192.3 13.45 226 -14.294  <.0001
 Grass E - S. bracteolactus       -209.4 14.22 226 -14.718  <.0001
 M. spinosum - S. bracteolactus    -17.0 11.71 226  -1.454  0.4671

P value adjustment: tukey method for comparing a family of 4 estimates 

# 2.3.2 - Parameter "w" ---- 
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")

 contrast                       estimate   SE  df t.ratio p.value
 Grass W - Grass E                  20.7 2.59 226   7.981  <.0001
 Grass W - M. spinosum             -31.7 7.75 226  -4.085  0.0004
 Grass W - S. bracteolactus        -27.0 8.07 226  -3.342  0.0053
 Grass E - M. spinosum             -52.4 6.62 226  -7.914  <.0001
 Grass E - S. bracteolactus        -47.7 6.84 226  -6.974  <.0001
 M. spinosum - S. bracteolactus      4.7 6.72 226   0.699  0.8974

P value adjustment: tukey method for comparing a family of 4 estimates 

# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time)/s))) + w, 
               start = c(A=15, m=30, s=-17, w=5),
               data = lfmc) # here, time is the only predictor variable 

b <- coef(nl.fe.0) # coefficients of the base model

# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) + w[leaf.type], 
              start = list(A=rep(b[1], 4), m=25, s=-10, w=rep(b[4], 4)),
              data = lfmc) # and fit is made using the base model's coefficients as start values

b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables

# 3.1.3 - Plot fixed effect ---- 
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group], 
               start = list(A=rep(c(b1[1],b1[2],b1[3],b1[4]),3), m=25, s=-10, w=rep(c(b1[7],b1[8],b1[9],b1[10]),3)),
               data = lfmc) # the model is fit using "b1" as the start values

IC_tab(nl.fe.0, nl.fe.1, nl.fe.2)

# Residual checking:
par(mfrow=c(1,1))
plot(nl.fe.2) # heterocedasticity in the residuals

hist(resid(nl.fe.2), main="", xlab="Residuals") # slightly skew distribution   

qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))

M2 <- nl.fe.2
summary(M2)

Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
A1    53.683      8.602   6.240 2.20e-09 ***
A2    54.282      8.636   6.285 1.72e-09 ***
A3    61.906      8.881   6.971 3.59e-11 ***
A4   111.561     13.291   8.394 5.65e-15 ***
A5   314.131     24.840  12.646  < 2e-16 ***
A6   333.831     26.640  12.531  < 2e-16 ***
A7    77.465     10.997   7.044 2.34e-11 ***
A8   298.161     24.917  11.966  < 2e-16 ***
A9   321.433     25.765  12.475  < 2e-16 ***
A10   87.553     11.744   7.455 2.01e-12 ***
A11  279.437     20.829  13.416  < 2e-16 ***
A12  292.443     24.178  12.095  < 2e-16 ***
m     30.413      2.896  10.504  < 2e-16 ***
s    -20.695      3.606  -5.739 3.11e-08 ***
w1    28.369      6.537   4.340 2.17e-05 ***
w2    27.480      6.548   4.197 3.93e-05 ***
w3    25.962      6.633   3.914 0.000121 ***
w4     0.890      8.173   0.109 0.913388    
w5    43.524     13.567   3.208 0.001535 ** 
w6    41.223     14.429   2.857 0.004686 ** 
w7     6.091      7.262   0.839 0.402466    
w8    26.610     13.604   1.956 0.051725 .  
w9    39.495     14.009   2.819 0.005252 ** 
w10    2.265      7.551   0.300 0.764455    
w11   66.661     11.478   5.808 2.19e-08 ***
w12   38.222     13.031   2.933 0.003708 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.62 on 221 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 7.269e-06

# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~1 | plot, data = lfmc)

# Residual checking:
plot(l.me) # heterocedasticity in the residuals

plot(ACF(l.me, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # temporal correlation is observed at lag 1

# 4.1.2 - Variance modelling ---- 
l.me.vm <- update(l.me, weights = varIdent(form=~ 1 | leaf.type))
    
# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled

AIC(l.me, l.me.vm) # modeling the variance improves the model fit 

# 4.1.3 - Temporal correlation modelling ----  
l.me.vm.tm <- update(l.me.vm, correlation=corARMA(p=1,q=1, form=~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # the temporal dependence is modeled

AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit

plot(l.me.vm.tm) 

plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) 

hist(resid(l.me.vm.tm, type="normalized"), main="", xlab="Residuals") 

qqnorm(resid(l.me.vm.tm, type="normalized"))
qqline(resid(l.me.vm.tm, type="normalized"))

# 4.1.4 - Evaluation of the fixed-effects ---- 
l.me.ml <- update(l.me.vm.tm, method ="ML") # the model is fitted by Maximum likelihood 
l.me.ml.1 <- update(l.me.ml, ~.-time:leaf.type) # the model is reduced removing the interaction "time x leaf type" 
IC_tab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important

M3 <- l.me.vm.tm 
summary(M3)

Linear mixed-effects model fit by REML
  Data: lfmc 

Random effects:
 Formula: ~1 | plot
        (Intercept) Residual
StdDev:    3.750488 7.842661

Correlation Structure: ARMA(1,1)
 Formula: ~1 | plot 
 Parameter estimate(s):
      Phi1     Theta1 
 0.9076001 -0.7567060 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W      M. spinosum S. bracteolactus          Grass E 
        1.000000         2.963571         3.123740         1.070897 
Fixed effects:  lfmc ~ time * leaf.type 

 Correlation: 
                               (Intr) time   lf.tGE lf.M.s lf.S.b tm:.GE
time                           -0.557                                   
leaf.typeGrass E               -0.709  0.395                            
leaf.typeM. spinosum           -0.407  0.227  0.654                     
leaf.typeS. bracteolactus      -0.394  0.219  0.653  0.666              
time:leaf.typeGrass E           0.374 -0.671 -0.592 -0.430 -0.429       
time:leaf.typeM. spinosum       0.176 -0.316 -0.299 -0.742 -0.409  0.546
time:leaf.typeS. bracteolactus  0.166 -0.299 -0.318 -0.441 -0.743  0.562
                               t:.M.s
time                                 
leaf.typeGrass E                     
leaf.typeM. spinosum                 
leaf.typeS. bracteolactus            
time:leaf.typeGrass E                
time:leaf.typeM. spinosum            
time:leaf.typeS. bracteolactus  0.567

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.3449313 -0.6438200 -0.1011963  0.6158869  3.3469739 

Number of Observations: 247
Number of Groups: 6 

M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)

Call:
lm(formula = lfmc ~ time * group, data = lfmc)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.673  -7.042   0.466   7.002  66.332 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          48.79275    6.44309   7.573 9.60e-13 ***
time                 -0.24478    0.13305  -1.840  0.06714 .  
groupGW plot 5        0.49817    9.11190   0.055  0.95645    
groupGW plot 6        6.14984    9.11190   0.675  0.50042    
groupGE plot 1       40.32675    9.49655   4.246 3.19e-05 ***
groupSM plot 1      212.68278    9.11190  23.341  < 2e-16 ***
groupSS plot 1      227.42468    9.11190  24.959  < 2e-16 ***
groupGE plot 2       14.39850    9.49655   1.516  0.13089    
groupSM plot 2      195.54780    9.11190  21.461  < 2e-16 ***
groupSS plot 2      217.12857    9.11190  23.829  < 2e-16 ***
groupGE plot 3       21.68513    9.49655   2.283  0.02334 *  
groupSM plot 3      189.79613    9.49655  19.986  < 2e-16 ***
groupSS plot 3      192.53490    9.49655  20.274  < 2e-16 ***
time:groupGW plot 5  -0.01905    0.18816  -0.101  0.91944    
time:groupGW plot 6  -0.10231    0.18816  -0.544  0.58718    
time:groupGE plot 1  -0.80129    0.19357  -4.139 4.94e-05 ***
time:groupSM plot 1  -2.36256    0.18816 -12.556  < 2e-16 ***
time:groupSS plot 1  -2.55767    0.18816 -13.593  < 2e-16 ***
time:groupGE plot 2  -0.43459    0.19357  -2.245  0.02574 *  
time:groupSM plot 2  -2.34722    0.18816 -12.474  < 2e-16 ***
time:groupSS plot 2  -2.45552    0.18816 -13.050  < 2e-16 ***
time:groupGE plot 3  -0.56657    0.19357  -2.927  0.00378 ** 
time:groupSM plot 3  -1.82099    0.19357  -9.407  < 2e-16 ***
time:groupSS plot 3  -2.16847    0.19357 -11.202  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.51 on 223 degrees of freedom
Multiple R-squared:  0.9586,    Adjusted R-squared:  0.9544 
F-statistic: 224.7 on 23 and 223 DF,  p-value: < 2.2e-16

# Residual checking:
par(mfrow=c(2,2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met) 

# 6 - NULL MODEL (M5) ######################################################################################

M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)

Call:
lm(formula = lfmc ~ 1, data = lfmc)

Residuals:
   Min     1Q Median     3Q    Max 
-88.75 -58.05 -31.45  52.66 231.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   95.645      4.919   19.45   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77.3 on 246 degrees of freedom

IC_tab(M2, M3, M4, M5)#, weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best

IC_tab(M1)

par(mfrow=c(2,2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")

# Nonlinear mixed-effects model (M1) 
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)  

boxplot(rM1 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)

plot(rM1 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)

qqnorm(rM1, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM1)

# Nonlinear fixed-effects model (M2)
par(mfrow=c(2,2), cex = 0.65) # (Figure 5)

rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)  

boxplot(rM2 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)

plot(rM2 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)

qqnorm(rM2, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM2)

par(mfrow=c(2,2), cex = 0.65)

# Linear mixed-effects model (M3) 
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)  

boxplot(rM3 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)

plot(rM3 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)

qqnorm(rM3, main="", ylab="Sample quantiles", xlab="Theoretical quantiles")
qqline(rM3)

par(mfrow=c(2,2), cex = 0.65)

# Clasical regression (M4) 
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)  

boxplot(rM4 ~ lfmc$leaf.type, ylab="Residuals", xlab="", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)

plot(rM4 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)

qqnorm(rM4, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM4)

—-